library(tidyverse) # Data manipulation + plotting libraries
library(sf) # Modern spatial vector data
library(nhdplusTools) # Access National Hydrography data (HUCs, flowlines, waterbodies)
library(elevatr) # Download elevation rasters
library(hillshader) # Generate hillshade from slope/aspect rasters
library(raster) # Masking/clipping rasters (used before terra conversion)
library(terra) # Terrain derivatives (slope/aspect), raster math
library(tidyterra) # ggplot2-friendly raster visualization
library(scales) # Rescaling + pretty axis/legend formatting
library(cowplot) # Compose final map + inset# load US state polygons
us_states <- maps::map_data("state")
# select only California
cali <- us_states[us_states$region == "california",]
# convert polygon coords to sf geometry
polygon <- sf::st_polygon(
list(as.matrix(cali[, c("long", "lat")]))
)
# create an sf object from polygon
cali_sf <- sf::st_sf(
id = "california",
geometry = sf::st_sfc(polygon),
crs = 4326
)# download HUC8 polygons
huc8 <- nhdplusTools::get_huc(
id = c("18020121", "18020122", "18020123"),
type = "huc08"
) %>%
st_transform(4326) # ensure consistent CRS
# dissolve multiple HUC8s into one combined polygon
feather_dissolve <- huc8 %>%
summarise(geometry = st_union(geometry)) %>%
st_as_sf()
# compute bounding box
feather_bbox <- st_bbox(feather_dissolve)
# convert bbox to polygon
feather_bbox_sf <- st_as_sfc(feather_bbox)# download NHD flowlines
nhd_flow <- nhdplusTools::get_nhdplus(
AOI = feather_dissolve,
realization = "flowline"
) %>%
st_transform(4326)
# Add custom linewidths by stream order
# cartographic hierarchy
nhd_flow <- merge(
nhd_flow,
data.frame(
streamorde = 1:6,
width = rev(c(1, .8, .6, .4, .2, .1))
)
)
# download lakes/ponds
nhd_wb <- nhdplusTools::get_waterbodies(
AOI = feather_dissolve
) %>%
st_transform(4326)# download DEM tiles
elevation_data <- elevatr::get_elev_raster(
locations = feather_dissolve,
z = 10, # zoom level
prj = st_crs(4326)$proj4string
)
# clip DEM to watershed
elevation_data <- raster::mask(
elevation_data,
feather_dissolve
)
# convert raster to terra object
r <- terra::rast(elevation_data)# hcl.colors() — grayscale palette for hillshade
pal_greys <- hcl.colors(1000, "Grays")
# Rescale hillshade values to palette indices
index <- hill %>%
mutate(
index_col = scales::rescale(shades, to = c(1, length(pal_greys)))
) %>%
mutate(index_col = round(index_col)) %>%
pull(index_col)
# Map grayscale values to palette
vector_cols <- pal_greys[index]
# draw rasters
map_hypso <- ggplot() +
tidyterra::geom_spatraster(
data = hill,
fill = vector_cols,
maxcell = Inf,
alpha = 1
) +
tidyterra::geom_spatraster(
data = r,
maxcell = Inf,
show.legend = FALSE
) +
tidyterra::scale_fill_hypso_tint_c(
limits = c(0, 2730),
palette = "wiki-2.0_hypso",
alpha = 0.6
) +
theme_minimal()# Add lakes/ponds
map_hypso_nhd <- map_hypso +
ggplot2::geom_sf(
data = nhd_wb[nhd_wb$ftype %in% "LakePond",],
fill = "darkslategrey",
color = NA
) +
# Add flowlines with hierarchical widths
ggplot2::geom_sf(
data = nhd_flow,
ggplot2::aes(linewidth = width),
color = "darkslategrey", show.legend = F
) +
scale_linewidth(range = c(0.05, .3))# Build inset map
inset <- ggplot() +
geom_sf(data = cali_sf, fill = "gray90") +
geom_sf(data = feather_dissolve, fill = "grey20") +
geom_sf(data = feather_bbox_sf, fill = NA, color = "black") +
theme_void()
# compose final map
final_map <- cowplot::ggdraw() +
cowplot::draw_plot(map_hypso_nhd) +
cowplot::draw_plot(inset, x = .73, y = .65, width = 0.25, height = 0.25)